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A novel scheme to solve the quantum eigenvalue problem through the imaginary-time Green lunction 
Monte Carlo method is presented. This method is applicable to the excited states as well as to the 
ground state of a generic system. We demonstrate the validity of the method with the numerical 
examples on three simple systems including a discretized sine-Gordon model. 
(to be published in Progress of Theoretical Physics 96 vol.5 (1996)) 
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The stochastic quantization has been successfully ap- 
plied in recent years to the various problems of physics 
ranging from quantum chromo-dynamics to the plasmon 
oscillation While it can handle very complex sys- 

tem which is usually beyond the approach either by the 
direct diagonalization or by the perturbation, it has been 
only applicable to the calculation of the ground state 
properties of the system. In this paper, we propose a 
new method based on the imaginary-time Green func- 
tion Monte Carlo approach, which is extended to the ex- 
cited states of Hamiltonian systems. We envision the 
application of the method both in non-relativistic quan- 
tum few-body problems [§]-f§ , and in lattice field theories 
0. The method is illustrated through the examples of 
a Morse oscillator, of a simple schematic deformed shell 
model with two interacting particles in two-dimension, 
and of sine-Gordon coupled oscillators. 

Let us suppose that we want to solve an eigenstate 
problem of a Hamiltonian H, namely 



H \ip a ) = E a \ip a ) 



(1) 



We assume that each eigenstate is normalized to unity. 
The evolution operator with imaginary time r acts on an 
arbitrary state fa as 



e- TH \<f> 1 )=Y,^ TEa \1> a )(1>a \fa}- 



(2) 



If we discretize the imaginary time r with the unit Ar, 
and call the state at the n-th step of the evolution 
Ifa)^ = exp — uAtH \fa), the evolution at each step 
Ar is described by 

|<M (n + l) =e -Ar ff|0i) («)_ (g) 

After sufficient number of steps, the state with the lowest 
energy is filtered out: 
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We consider the second sequence of states V2 
satisfy the evolution equation 
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The states <j)^ and fa^ 1 ' at adjacent time steps satisfy the 
orthogonality ^{fa\ fa) = 0. Expanding eq. (||) in 



the first order of Ar, we have 



v(n+l) 
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whose Ar-dependent term in the first line vanishes at 
the large time step because of eq. It is now clear 

that the repeated operation of eq. (||) filters out the first 
excited state, namely: 
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It is important to notice that the orthogonalization at 
each time step, eq. (||) is essential in actual calcula- 
tion, because the single Schmidt subtraction either at the 
starting or at the end of the repeated iteration is imprac- 
tical with finite accuracy numerics. One can generalize 
eqs. (||) and (|J) to obtain the evolution equation for the 
a-th energy eigenstate (a — 1, 2, 3, ...) as 



\fa 



- ArH \fa) {n) (8) 
{n) {fa\e- ATH \fa) in) 
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Then, with the step-wise asymptotic orthogonal condi- 
tions W (0^1 <j) a )^ n+1 ^ — (n — oo) for different states 
a and (3 , one can recursively filter out the a-th energy 
eigenstate as 



1 



\lpa) (lpa\ <t>a) (n — ► OCi). (9) 



We work in configuration space, and adopt the usual 
notation (f>a (q) = (q\ <t>a)^ where q signifies the co- 
ordinate vector of the problem. The straightforward 
evaluation of integrals in the evolution equation is im- 
practical except for the case of very low dimension. In- 
stead, we consider, for given n and a, a set of vectors 
{q^}(i = 1...M) distributed according to the probabil- 
ity distribution 



P(q) = 



dq ^(q) 



(10) 



and try to represent the wave function by the method of 
importance sampling. The sign information of <p^{q), 
which is missing in eq. ( |l0| ) , can be recovered by defining 
a real number 



» 



sgn 



(11) 



dq 



for each point q^ . One can evaluate any inte- 
gral involving the wave function (jy^iq) using the set 

Wai\ s ai} (* = ^-■■■M). The most relevant example is 
the coordinate representation of the evolution equation 
eq. (ph, which now reads 



^ +1) (q) = / dq'K Ar (q,q")^\q) 



(12) 
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where K Ar (q, q') = (g| exp(— AtH) \q') is the step prop- 
agator in g-representation. The overlap matrix A^j 7 is 
given by 

A^ = W^le-Arff^^Cn) (13) 

Another overlap matrix A,g 7 = ^(4>(j\ </> 7 /™ +1 ^ is calcu- 
lated with the recursion relation 



7-1 



a=l 



(14) 



Since we always take care of the normalization of the 
wave functions explicitly, the overall factor in the inte- 
grals of eq. (|l^) is irrelevant, and one can discard the un- 
seemly constant (integral and sum) in the definition eq. 
(|ll|) to regard as simple sign, = ±1. The stan- 
dard Metropolis algorithm M can be applied on eq. ( p^ ) 
to obtain a new set of numbers which 
samples the 4^a +1 \q) from the set {q^ , s£™ '} by a ran- 
dom walk. At each step, one can evaluate the energy 
integral with the sampling points and associated sign as 



_ Jdq^ a n> (q)H^ a n+L >(q) 
: Jdq^\q)^ +1) (q) 

EHK AT (q£\qff) S ^s% 



(15) 
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where HK Ar (q,q') = (q\ Hexp(-ArH) \q' ). This se- 
quence of integrals, after the convergence, which we as- 
sume to occur before the no-th step, should yield the true 
eigenvalue 
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After the n -th step, all N s * M sets {q^> ,s^' } (i = 
l...M,n — nQ...n + N s ) can be thought of as sampling 
the true eigenstate wave function that can be used to 
calculate any observable. It should be noted that the eq. 
( |l2| ) is exact only in the first line, and the second line 
has the sampling error in the estimation of the overlaps, 
and includes the contamination by the lower eigenstates. 
However, because of the law of large numbers, the fluc- 
tuating error at each step should cancel each other in a 
long run, as long as the sequence of the state <p& con- 
verges to a fixed limit. Therefore, in large N s limit, E a 
of eq. ( |i~6| ) does converge toward the true eigenvalue for 
any excited states, albeit more slowly for higher states. 

Now we turn to the numerical examples. Here, we 
can only outline the results deferring the full detail to 
a forthcoming publication. First, we look at the one- 
dimensional motion of a particle in a Morse potential to 
see the workings of the method in a simplest setting. We 
take the Hamiltonian H 

TABLE I. The energy eigenvalue of the first five states of 
the Morse system. The raw [a] is the results with At = 0.5, 
M = 200 and N = 400, while [b] is with At = 0.2, M = 1000 
and N = 80. 



a 


1 


2 


3 


4 


5 


E [a] 


-7.00 


-5.29 


-3.90 


-2.86 


-1.84 


E [b] 


-7.01 


-5.27 


-3.79 


-2.56 


-1.53 


(Exact) 


(-7.031) 


(-5.281) 


(-3.781) 


(-2.531) 


(-1.531) 
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H = --—+8e~ q/2 -16e~ q . (17) 
2 oq l 

In Table 1, the first five eigenvalues evaluated with two 
sets of sampling parameters (M = 200, N s = 400, At = 
0.5) and (M = 1000, N s = 80, At = 0.2) are tabulated. 
Statistical errors are estimated to be 0.2% at most for 
both cases. One can see, in the table, that the results 
do converge toward the analytical value E a = — (17 — 
2a) 2 /32 with larger M and smaller At. 
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FIG. 1. The profile of the first five wave function calcu- 
lated with Monte Carlo technique. 

In Fig. 1, we plot the wave functions, which are ob- 
tained as suitably normalized probability distributions 
of sampling points taken separately for s^i = 1 and 

(n) 

s ai — — 1. Their satisfactory quality is visually evident. 

We look at another numerical example of two iden- 
tical particles in a common two-dimensional asymmet- 
ric harmonic potential, mutually interacting through a 
Gaussian force. This can be regarded as a simplified, but 
non-trivial prototype of the deformed shell model [fTo) . 
The model Hamiltonian is 

f (gi - x 2 f + {yi- V2? \ 
+,oexp| - j. 

The coordinate vector q = (g*(l) , g*(2) ) = (xi, y\, x 2l y 2 ) is 
four dimensional. We assume the particles to be fermions 
of 1/2 spin. The spin degree of freedom and the identity 
of the particle pose us the question of exchange symme- 
try. We consider a generic S z = state which we write 
as 

<M#), Q^)uid 2 + cf> du (qXl), q(2))d lU2 (19) 



where u and d represent the S z = 1/2 and Sz = —1/2 
states of individual particles respectively. From eq. (|l^) , 
we can construct the totally antisymmetric state 

{<f> u d(<ffl, q(2)) - cf> du (q[2),q{l))} "l* (20) 
+ {^ du (q\l), q(2)) - 4> ud (q\2),q\l))} d x u 2 

which contains all the physical states. Let us assume 
that the wave functions <j) ud and 4>d u are sampled respec- 
tively by the sets {q u d(l)i, $ud®i, s u di} (i = 1...M) and 
{qdu($i,qdu(2)i,s du i} (i = 1...M). Clearly, one can sam- 
ple each term of eq. (20) by enlarging the sets to include 
2M elements with appropriate signs: Namely, for the 
first term, {q u dQ)i, q-udtyi, Sudi} © {<fd«(2)i, -Srf™} 
(i = 1...M), and for the second, {q du (i)i, qdu®i, Sdui} ® 
{q u d($i,qud(l)i,-Sudi} (i = l—M). The operation of 
enlarging the sampling points for antisymmetrization 
should be carried out at each step of the evolution in eq. 
(|l2|). Otherwise, the system would pickup the small con- 
tamination of other symmetries, and quickly evolve into 
the lowest energy state of any symmetry, namely the to- 
tally symmetric state. This scheme for the treatment of 
the exchange symmetry is generalizable, in principle, to 
the system with arbitrary number of particles. 

The results of the calculation for the Hamiltonian eq. 
( |l8| ) is summarized in Table 2. The values for the inter- 
action parameters are taken to be vq = 1 and a = 0.5. 
The sampling number 2M is set to 1600. The time step 
is chosen to be At = 0.5. The sampling point from N s — 
40 iterations (after the initial iterations of several hun- 
dred) are included. The first line is the energy eigen- 
value of first five states. The second line lists the value 
of 45(5 + 1) which should serve as a barometer for the 
quality of wave functions with respect to the separate 
spatial and spin symmetries: It should be for spatially 
symmetric 5 = states, and 8 for spatially antisymmetric 
5=1 states. It appears that the third state still contains 
the contamination by components of wrong symmetry by 
almost 10%. This suggests that the choice of either At 
or 2M is still insufficient to obtain the full convergence. 
But the overall feature of the eigenstates are reasonably 
well captured even at this level. The one-body density 
profiles calculated from the sampled points are depicted 
in Fig. 2. Because of the repulsive two body force, the 
spatially symmetric states are pushed up in the energy 
compared to the antisymmetric companion states. 

TABLE II. The energy eigenvalue and the total spin of 
the first five states of the Gaussian coupled two particles in 
two-dimensional oscillators. The truncation parameters are 
At = 0.5 and 2M = 1600. 
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2 


3 


4 
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E 


2.44 


3.34 


3.43 


3.56 


3.69 


4S' 2 


0.11 


7.56 


0.73 


7.27 


0.28 
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FIG. 2. The one body density profiles of the first five 
eigenstates of the Hamiltonian eq. (18). 

Finally, we look at a linear array of oscillators whose 
Hamiltonian is given by 



H = 



N 

E 

fc=i 



1 

2a 5 ' 



(21) 



9 N 



(27) 



are calculated with the iterated Monte Carlo sampling as 
before. The only complication is that we need to include 
all ./V-cyclic permutations in the set of sampling points 
at each step of the iteration to ensure the projection, 
eq. (^6|): Namely, if ^ nf ,a is sampled by a set of M vec- 
tors {4>u, 4>2i—<pNii s i} {i — l.-.M), ^ nf ,a can be sampled 
by NM vectors {<j>u, (j) 2ll <j>m, Si} © {<f>2i, <j>u, —,<f>u + 
nf2z,Si}®...®{4>m, ^u+nf^-, <j>N-\%+nfjf, s t } (i = 
1...M). From the eigenvalues E nf ^ a , we can extract the 
mass of bosons as 



rribi = E . 2 — E i, mvi = E , 3 — Eq.i, 



etc. 



(28) 



where m^i is the mass of the elementary boson, m^, the 
second boson, etc.. The mass of the fermion is given by 



Mj = E\fi — Eq,o- 



(29) 



where </>& and life = <j>k are the amplitude and the conju- 
gate momentum of the fc-th oscillator, and a is the com- 
mon distance between the neighboring oscillators. The 
coupling between adjacent oscillator is given by a trigono- 
metric form 



Try 



(22) 



This is a discretized version of the one-dimensional sine- 
Gordon field. The eigen spectra with the periodic bound- 
ary condition 



4>N+1 = 01 



(23) 



give the bosonic excitation modes. There are also soliton 
modes which are identified as the fermionic excitations. 
They are given by the twisted boundary condition 



J N+1 



2n 
J 



(24) 



where an integer nt is the number of kinks, or fermions. 
Direct extraction of excitation spectra from the Hamil- 
tonian, eq. ( pl| ) is difficult, because of the existence of 
translational zero mode of the system. We define an op- 
erator T which shifts the whole system by a as 



T"a*(01, 4>2i ■■; <f>N-l, 4>n) 
*(</>2,<fe, 



(25) 



'■>N, 



27T 

ns T ] 



where VP is an arbitrary wave function. Clearly, T com- 
mutes with the Hamiltonian H. Therefore, from an ar- 
bitrary eigenstates ^ ni ,a of H, one can project out the 
translationally invariant component as 



n f , 



{l+T a +T 2 a 



T 



N-l 



nj, 



(26) 



The projected eigenstates, which are the solutions of the 
eigenvalue problem 
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FIG. 3. (a) The soliton mass Mf as the function of the 
coupling j3 and (b) the bosonic masses rribi and rrihi as the 
function of soliton mass Mf in one dimensional sine-Gordon 
model. 

In Fig. 3 (a), we plot the fermion mass Mf calculated 
with the Green function Monte Carlo procedure as the 
function of the coupling constant (3. In Fig. 3 (b), the 
boson masses rribi and mb2 are plotted as the function of 
the fermion mass Mt, The time step is taken to be At = 
0.5. The sampling number NM and the iteration number 
N s are chosen to be 360iV and 400 respectively. Results 



4 



for N — 4,6,8 are shown in the figures. In Fig. 3 (a), 
results for different N are practically indistinguishable. 
For the continuum limit N — > oo, the mass spectra of 
the sine-Gordon model has been calculated within semi- 
classical approximation using the trace formula . This 
analytical result is shown in the figures as the line. The 
agreement with the numerical results seems to support 
the assertion by the authors of rcf. p| that their results 
are exact in full quantum sense. 

In above examples, the computational time required 
for the convergent calculation still tends to exceed the 
time required for the conventional direct diagonalization 
approach. However, the expected increase in time and 
storage for lager degrees of freedom is more modest in the 
Monte Carlo calculation than in conventional approaches. 
Also, the Monte Carlo methods have inherent affinity to 
the computational parallelism. We therefore believe that 
the results obtained here are sufficiently encouraging for 
the prospect of tackling more realistic problems both in 
many-body quantum theories and field theories. 

We express our gratitude to Drs. O. Morimatsu, M. 
Takizawa, and T. Kawai for the helpful discussions. 
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